%==========================================================================
%=  This file is part of the Sound Restoration Project
%=  (c) Copyright Industrial Mathematics Institute
%=                University of South Carolina, Department of Mathematics
%=  ALL RIGHTS RESERVED
%=
%=  Author: Borislav Karaivanov
%==========================================================================

%==========================================================================
% List of the files on which this procedure depends:
%
% createMissingDirs.m
% setLogWritingFunctHandle.m
% visualizeOverlapGuess.m
% writeToLogFile.m
%
%==========================================================================

%==========================================================================
% The function "computeOverlapGuess" takes a collection of sound signals of
% the same length and computes the most popular amount of overlap among a
% specified number of pairs of consecutive signals chosen randomly among
% all available pairs of this kind.
% INPUT: "allSignalsArr" is a 2D array whose columns hold the individual
% signals whose overlaps are to be computed.
% "numToUse" is an optional, positive integer specifying the number of
% pairs of consecutive signals whose overlaps are to be computed. These
% pairs are to be chosen randomly among all available pairs of consecutive
% signals. Note that the actual number of used pairs may be slightly
% different due to possible repetitions which are removed.
% "locMinNbhdRad" is an optional, positive integer specifying the radius of
% the neighborhood (an interval, in fact, in 1D) around a point to be
% checked in order to claim a local minimum at that point. More precisely,
% we consider having a local minimum of the relative differences at a fixed
% overlap value if for all overlaps in the neighborhood of that value the
% relative differences are not smaller.
% "writeToLogFunctHandle" is an optional function handle which writes to a
% fixed log text file if keeping a log is requested, or does nothing
% otherwise.
% OUTPUT: "mostPopularOverlap" returns the most popular among the overlaps
% of the chosen pairs of consecutive signals. It is measured in number of
% samples.
% "percentPopularity" returns the percentage of chosen pairs whose overlap
% equals the most popular overlap.
%==========================================================================
function [mostPopularOverlap, percentPopularity] = computeOverlapGuess( ...
    allSignalsArr, numToUse, locMinNbhdRad, writeToLogFunctHandle)

% Get the number of signals and their common length.
[signalLength, numSignals] = size(allSignalsArr);

% Define default values for the optional input arguments.
if (nargin < 2)
    numToUse = 400;
end
if (nargin < 3)
    locMinNbhdRad = 5;
end
if (nargin < 4)
    writeToLogFunctHandle = @(x) '';
end

% Make sure that the specified number of signals to be used does not exceed
% the number of all signals.
numToUse = max(0, min(numToUse, numSignals - 1));
if (numToUse == 0)
    mostPopularOverlap = 0;
    percentPopularity = 100;
    return;
end

%=========================================================================%
%========== Here is the essential portion of this function. ==============%
%=========================================================================%
% Randomly pick the specified number of pairs of signals whose overlaps are
% to be computed, and remove any duplicates.
toUseIndArr = unique(randi([1 (numSignals - 1)], numToUse, 1));
% Get the number of pairs remaining after duplicates are removed.
numToUse = numel(toUseIndArr);

% Extract the signals that will participate in the computation of overlaps.
% Note that those portions are stored transposed because it is faster to
% reduce all signals by one sample at once if those samples are stored in a
% column.
firstSignalArr = allSignalsArr(:, toUseIndArr).';
secondSignalArr = allSignalsArr(:, toUseIndArr + 1).';

% Allocate memory for the average absolute differences corresponding to all
% overlaps and all chosen pairs of signals.
aveDiffArr = zeros(numToUse, signalLength);
aveSumArr = zeros(size(aveDiffArr));

% Set the number of pairs of images to be processed at once. The number 100
% is not picked arbitrary. It is the round number for which the code seems
% to perform best on the machine where tested. It hits the balance between
% the number of iterations and the size of the arrays used on each
% iteration.
numPairsAtOnce = 100;

% Compute the average absolute differences corresponding to all allowed
% overlaps and all chosen pairs of signals.
for pairInd = 1:numPairsAtOnce:numToUse
    
    % Get the indices of the pairs of consecutive signals in the current
    % batch to be processed.
    pairIndArr = (pairInd:min(numToUse, (pairInd + numPairsAtOnce - 1)));
    
    % Extract the portions of the signals in the current batch of pairs.
    currFirstSignalArr = firstSignalArr(pairIndArr, :);
    currSecondSignalArr = secondSignalArr(pairIndArr, :);
    
    for ind = signalLength:(-1):1
        
        % Compute the average absolute differences and sums corresponding
        % to the current batch of pairs of signals and the current overlap.
        aveDiffArr(pairIndArr, ind) = mean(abs(currFirstSignalArr - currSecondSignalArr), 2);
        aveSumArr(pairIndArr, ind) = mean(abs(currFirstSignalArr) + abs(currSecondSignalArr), 2);
        
        % Reduce the signals in the current batch of pairs by one sample.
        currFirstSignalArr = currFirstSignalArr(:, 2:end);
        currSecondSignalArr = currSecondSignalArr(:, 1:(end - 1));
    end
end

% Compute the relative average absolute differences. Below, for the sake of
% brevity, the relative average absolute difference will be called relative
% difference. It is important to use relative differences because otherwise
% a difference could be small just because two very low (in magnitude)
% portions of the signals are compared (there all samples have values close
% to zero and therefore the differences of the sample values can not be
% large). When relative differences are used we divide the usual difference
% by the sum of the L_1 norms of the two compared portions. Thus if the two
% portions are very low in magnitude we will be dividing by a small number
% and the relative difference will end up being larger as we would like it
% to be.
relDiffArr = aveDiffArr./aveSumArr;

% Find all local minima of the relative differences for each chosen pair of
% consecutive signals.
minRelDiffArr = ordfilt2(relDiffArr, 1, ones(1, 2*locMinNbhdRad + 1), 'symmetric');
isLocalMinArr = (minRelDiffArr == relDiffArr);

% For every allowed overlap compute the number of pairs having a local
% minimum of the relative average difference at that overlap.
countArr = sum(isLocalMinArr, 1);

% Get the maximal count and the position where it is attained. Note that 
% the first "locMinNbhdRad" and the last "locMinNbhdRad" possible overlaps
% are excluded from the considerations for the maximal count. This is
% motivated by the following two arguments: when the amount of overlap is
% too small, a large count might be indicating lack of sufficient
% discriminating features in the overlapping portions rather than
% well-matching features; when the amount of overlap is too close to the
% image height, a large count might be indicating resemblance between two
% consecutive frames rather than well-matching portions of two scans
% corresponding to one and the same frame.
firstInd = locMinNbhdRad;
lastInd = signalLength - locMinNbhdRad;
[maxCount, maxCountInd] = max(countArr(firstInd:lastInd));

% Compute the most popular overlap and its popularity (i.e., the percentage
% of pairs whose relative differences have a local minimum at the most
% popular overlap).
mostPopularOverlap = maxCountInd + firstInd - 1;
percentPopularity = 100*maxCount/numToUse;

% Get the "second largest" count, i.e., the largest count remaining when
% excluding a neighborhood of the first largest count.
[secondMaxCount, secondMaxCountInd] = ...
    max(countArr([(firstInd:(firstInd + maxCountInd - 1 - locMinNbhdRad - 1)), ...
    ((firstInd + maxCountInd - 1 + locMinNbhdRad + 1):lastInd)]));

% Compute the second most popular overlap and its popularity (i.e., the
% percentage of pairs whose relative differences have a local minimum at
% the second most popular overlap).
if (secondMaxCountInd <= maxCountInd - locMinNbhdRad - 1)
    secondMostPopularOverlap = secondMaxCountInd + firstInd - 1;
else
    secondMostPopularOverlap = secondMaxCountInd + firstInd - 1 + 2*locMinNbhdRad + 1;
end
percentPopularityOfSecond = 100*secondMaxCount/numToUse;
%=========================================================================%

% % Visualize and save snapshots, if desired.
% visualizeOverlapGuess(countArr, relDiffArr, isLocalMinArr, numToUse, ...
%     locMinNbhdRad, mostPopularOverlap, percentPopularity, secondMostPopularOverlap, ...
%     percentPopularityOfSecond);

% Write info to a log file.
writeToLogFunctHandle(['number of signals used for overlap guess: ' num2str(numToUse)]);
writeToLogFunctHandle(['local min neighborhood radius for overlap guess: ' num2str(locMinNbhdRad)]);
writeToLogFunctHandle(['best overlap guess: ' num2str(mostPopularOverlap)]);
writeToLogFunctHandle(['percent popularity of best overlap guess: ' num2str(percentPopularity, '%.2f')]);
writeToLogFunctHandle(['second best overlap guess: ' num2str(secondMostPopularOverlap)]);
writeToLogFunctHandle(['percent popularity of second best overlap guess: ' num2str(percentPopularityOfSecond, '%.2f')]);
writeToLogFunctHandle(['percent second best overlap guess is weaker: ' ...
    num2str(100*(percentPopularity - percentPopularityOfSecond)/percentPopularity, '%.2f')]);

return;
% end of the function "computeOverlapGuess"
